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We investigate solitons and nonlinear Bloch waves in Bose-Einstein condensates trapped in op- 
tical lattices. By introducing specially designed localized profiles of the spatial modulation of the 
attractive nonlinearity, we construct an infinite number of exact soliton solutions in terms of the 
Mathieu and elliptic functions, with the chemical potential belonging to the semi-infinite bandgap 
of the optical-lattice-induced spectrum. Starting from the exact solutions, we employ the relaxation 
method to construct generic families of soliton solutions in a numerical form. The stability of the 
solitons is investigated through the computation of the eigenvalues for small perturbations, and also 
by direct simulations. Finally, we demonstrate a virtually exact (in the numerical sense) composition 
relation between nonlinear Bloch waves and solitons. 
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I. INTRODUCTION 

Analogies between the electron dynamics in perfect 
crystals and light propagation in periodic optical media 
suggest a variety of physical phenomena and related ap- 
plications. Bose-Einstein condensates (BECs) in optical 
lattices (OLs) not only represent an ideal tool for inves- 
tigating fundamental effects, such as the Landau-Zener 
tunneling, Josephson oscillations, dynamical instabilities, 
and quantum phase transitions between the superfluidity 
and the Mott insulator, but also offer versatile setups for 
the potential implementation of quantum computation 
schemes P, @j. 

The mean-field description of the BEC dynamics at 
zero temperature is based on the Gross-Pitaevskii equa- 
tion (GPE), that is, the nonlinear Schrodinger equation 
(NLSE) with a potential term, which is a ubiquitous 
model with important realizations in other fields - first of 
all, nonlinear optics [3| . Many experimental and theoret- 
ical works P-flq (see also reviews [H,l2,[l7|,[3) have been 
dealing with matter- wave and optical solitons in OLs. 
Usually, these solitons are found in a numerical form, 
with their chemical potential falling into bandgaps of 
the spectrum induced by the OL potential, in the frame- 
work of the corresponding linear Schrodinger equation. A 
specific dynamical phenomenon, which is relevant to the 
present work, is the composition relation between nonlin- 
ear Bloch waves (NBWs) and fundamental gap solitons, 
whose main peaks are confined to a single OL cell [l9[ . 

Current experiments with BECs wield a high degree 
of control over key parameters of the system. By means 
of the Feshbach resonance-technique, driven by magnetic 
or optical fields [2(J Hl| , one can adjust almost at will 
the strength and sign of the inter-atomic interaction. On 
the other hand, available fabrication technologies allow 
a modulation of nonlinearity in nonlinear optics. There- 
fore, there has been increased interest in the study of the 
nonlinear dynamics under spatially modulated nonlinear- 



ities, in optics and BEC alike, see original works |22H35l | 
and book [37] . In such settings, the nonlinear dynam- 



ics exhibits novel features, such as the "anti-Vakhitov- 
Kolokolov" criterion which controls the stability of gap 
solitons in media combining a spatially periodic nonlin- 
earity and the OL potential 30] . 

Exact solutions for matter-wave solitons in BECs with 
OL potentials are important not only because of their 
simplicity and the connection to physical bound states, 
but also since they can be used to test various approx- 
imate methods, and may also find applications in other 
fields. The objective of the present work is to construct 
one-dimensional soliton solutions in physically relevant 
situations combining the OL potential and a spatially 
modulated attractive nonlinearity. In addition to produc- 
ing exact soliton solutions in specially devised versions 
of such systems and exhibiting their relation to NBWs, 
we also find generic numerical solutions, by means the 
relaxation method, and investigate their stability. The 
results may be also be directly applied to nonlinear op- 
tical media with embedded periodic gratings, which play 
the same role in photonics as the OLs in BEC. 



II. THE MODEL AND ITS REDUCTION 

We consider a condensate of atoms trapped by a com- 
bination of a tight cigar-shaped magnetic trap and an 
OL potential acting in the longitudinal direction. If 
the transverse dimensions are comparable to the healing 
length, and the longitudinal dimension is much longer 
than the transverse ones, the setting is effectively one- 
dimensional, obeying by the respective version of the 
GPE (see, e.g., Refs. @): 



-yj xx + [2V cos(2x)+g(x)\^\ 2 },p, (1) 



where vj(x,t) is the macroscopic wave function of the 
condensate. Here, time t, spatial coordinates x, and 



2 



the strength of the OL potential, Vq, are normalized, 
respectively, by K/E r , k, and E r /A, with the recoil en- 
ergy E r = h 2 k 2 /2m, wave number of the optical lattice 
k, and atomic mass to. The nonlincarity coefficient is 
g = Amu! r a s /hk 2 , where u r is the transverse harmonic 
frequency and a s is the s-wave scattering length of inter- 
atomic collisions. By means of the Feshbach-resonance 
technique controlled by properly designed configurations 
of external fields, a s may be subject to a spatial modu- 
lation, hence the corresponding nonlinearity coefficient, 
g(x), may be a function of x. In this paper, we focus on 
the attractive nonlinearity. namely, g(x) < 0, rather than 
more general situations with the sign-changing g(x), such 
as those considered in some other works (see, in partic- 
ular, Refs. p9l. [3fl l35l] ) . It is relevant to mention that 
the cubic nonlinearity in Eq. (TTJ) is valid if the density 
is small enough; otherwise, the reduction of the dimen- 
sion in the GPE from 3 to 1 leads to a nonpolynomial 
nonlinearity [36| . 

Stationary soliton solutions to Eq. (TT]) are searched as 
rp(t,x) = (j>(x) exp(— i/Mt), where chemical potential p is 
normalized by the recoil energy, and real function <fi(x) 
obeys the following stationary NLSE, 



flip = -<p xx + [2V cos(2a;) + g(x)(f) 2 ](b, 



(2) 



with boundary conditions <p(x — > ±00) = 0. Up to the 
rescaling, the number of atoms and energy of the local- 
ized state are 



N = 



dx, 



(3) 



E„ 



fiN 



dx 



+ 2V cos(2a;)0 2 + ^-(f> 



dx 



(4) 



Following the scheme proposed in Refs. [26ll27|. exact 
soliton solutions can be constructed by casting Eq. ([2]) 
into the from of a solvable stationary NLSE in the free 
space, 



EU = -U X x+g U J 



(5) 



where E and go are constants. This reduction may be 
implemented by employing the transformation, 



It is commonly known that Eq. ([5]) possesses exact so- 
lutions in terms of the Jacobi's elliptic functions. There- 
fore, exact soliton solutions to Eq. ([2j can be constructed 
as long as exact solutions of Eq. © are known. In fact, 



P : 



aip 2 



(9) 



solves the Ermakov-Pinney equation, where a, (3 and 
7 are real constants satisfying E — (cry — /3 2 ), ip\ = 
MathieuC(/^, Vq, x) and ^2 = MathieuS(/i, Vb, x) are two 
linearly independent Mathieu functions that satisfy the 
Mathieu equation [H, (p xx + [m — 2Vb cos(2x)]<y5 = 0. 
For the soliton solutions to be physical, from Eq. (J7J) it 
follows that p(x) must not change its sign at any point 
(a sign-definite function), otherwise the local nonlinearity 
would diverge at points of p = 0. Therefore, parameters 
a, (3, 7 in Eq. ^ should be chosen so as to secure this 
condition. 



III. EXACT SOLITON SOLUTIONS WITH THE 
ATTRACTIVE NONLINEARITY 

For the attractive nonlinearity, go < 0, a relevant exact 
nontrivial solution to Eq. ([5J) is 



U(X) = ^(E~X 2 )/g cn(XX - X , to), (10) 

where A and Xq are two arbitrary constants, E satisfies 
— A 2 < E < A 2 , and cn is the Jaco bi's elliptic function 
with module to = ^(A 2 - E)/2X 2 . When \E\ < A 2 , 
Eq. (1101) gives a periodic function of X, with the min- 
imum period 4K(m)/\X\, where K(m) is the complete 
elliptic integral of the first kind. Since p(x) 7^ 0, the 
boundary condition <p(x — ¥ ±00) = is satisfied when 
U(X(x — > ±00)) = 0. According to the periodicity 
of Jacobi cn function, X[X(x —> +00) — X(x — > —00)] 
should be 2nK(m), with integer n. When E = —A 2 , 
solution (jlOl) goes over into the well-known elementary 
one, U(X) -- 



y / 2E/g sech {\/—EX — Xq\, where the 
boundary condition, <p(x —> ±00) = 0, may be satisfied 
if X(x — > ±00) = 00. 



A. The case of E = 

For E = 0, Eq. (jSJ is linear, and its solution can be a 
linear combination of the Mathieu functions, 



x) = p(x)U[X(x)], X(x)= p{s)~ 2 ds, (6) p = Cl MathieuC(/i,F ,a;) + c 2 MathieuS(/x, Vb, z), (H) 



and requiring 

g{x) = go P - & (x), (7) 
where p obeys the Ermakov-Pinney equation [26-28, 38], 
p xx + [p- 2V cos(2x)]p = Ep~ 3 . (8) 



where the constants c\ and C2 should be chosen so as to 
make p (x) sign-definite. 

We begin by constructing exact symmetric and anti- 
symmetric soliton solutions for Eq. ([2]) , where the spatial 
modulation of the nonlinearity should be represented by 
an even function p{x). Without the loss of generality, we 
then set c\ = 1 and C2 = in Eq. (jlip . hence p is an even 
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FIG. 1: (a) The cutoff chemical potential versus the strength 
of the OL potential (solid line). Exact soliton solutions exist 
below the solid line. The circle designates p co = —4.2805 at 
Vo = 4. (b) The spatially-modulated nonlinearity coefficient, 
as given by Eq. (fT5)) . with p = —10, Vo = 4, and go = —1. 
(c) The same as (b), except p = —4.2807. 



FIG. 2: (Color online) Exact symmetric solitons for (a) n = 1 
and (e) n = 3, and an exact antisymmetric soliton for (c) 
n — 2, where the corresponding modulation function g(x) is 
taken as per Fig. 1(b). Panels (b), (d) and (f) are the same 
as (a), (c), and (e), respectively, expect the corresponding 
modulation function is taken as in Fig. 1(c). Solid circles in 
the bottom of each column show the OL potential. 



function of x. Since p{x) should also be a sign-definite 
function, chemical potential p cannot be arbitrary for 
fixed strength Vq of the OL potential. It can then be 
shown that, for given Vo, there is a cutoff value of the 
chemical potential, namely, 



p co = MathieuA(0, Vo), 



(12) 



below which p is sign-definite. Here p co is an even func- 
tion of Vo, representing the first characteristic value of the 
MathieuC function, so that MathieuC is a 27r-periodic 
function of x. The Taylor expansion for small Vq is p co 
= -(l/2)Vo 2 + (7/128)Vb 4 + 0(Vb 6 ), with p co = at 
Vb = 0. The cutoff chemical potential versus Vq is shown 
in Fig. 1. 

Interestingly, we find that p co is exactly the minimum 
energy eigenvalue in the first Bloch band of the corre- 
sponding liner Schrodinger equation with periodic po- 
tential 2Vb cos(2a;). Thus, these exact soliton solutions 
of Eq. {2} exist in the semi-infinite bandgap. 

Now we investigate the properties of the nonlinearity- 
modulation pattern and respective solitons. When p is 
much smaller than the cutoff value p co , p(x) increases 
monotonically and quickly approaches infinity. There- 
fore, the modulation function, 



.9o 



MathieuC (p, Vb, x) 



(13) 



is localized in a very narrow single region [Fig. 1(b)]. 
Also, from Eqs. (J7|) and (jSJ it follows that the smaller the 
chemical potential, the narrower the localization region. 
On the contrary, when p approaches p co , p oscillates and 
slowly approaches infinity, so that the region of the local- 
ization of g(x) is relatively wide, featuring several layers 
[Fig. 1(c)], and the more closely the chemical potential 



approaches the cutoff value, the wider the localization 
region of the nonlinearity coefficient. 

Since the even and sign-definite p(x) approaches in- 
finity at \x\ — > oo, it is clear that X(x), defined in 
Eq. (|6]), is a monotonic non-decreasing odd function of 
x, which has upper and lower limits. Therefore, to let 
the exact soliton solutions meet the boundary condition 
<fi(x — > ±oo ) = 0, constant A in Eq. ([TO]) must be chosen 
so as to satisfy condition XX (x — > +oo) = nK{\2/2), 
where n = 1,2, 3, At the same time, constant Xq 
should be chosen as Xq — for even integer n, and 
X = K(l/V2) for odd inte ger n. Thus, exact soliton 
solutions to Eq. ([2]). with the modulation pattern taken 
as per Eq. (fl"3"]). are 



(t>ni.x) = 



x cn 



^9o"A(+oo) 
/ nK(l/V2) 
{ X(+oo) 



MathieuC(^i, Vq, x) 



X,l/V2\ , 



(14) 



for n — 1, 3, 5, ... , while for n = 2, 4, 6, ... the exact solu- 
tions are 



<l>n(x) 



x cn 



/-goX(+oo) 

nK (1/V2) 



MathieuC (p, Vq, x) 



X - K [ — - ] > — — 

A(+oo) I 2 I 2 



(15) 



where we define X(x) = MathieuC(^, Vq, s)~ 2 ds. 

It follows from Eqs. ©, (JTSJ), 03]), and JTS} that, 
once the chemical potential (p < p co ), constant an d 
the strength of the OL potential, Vq, are fixed, there ex- 
ists an infinite number of exact solitons sharing the same 
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FIG. 3: (Color online) Exact symmetric and anti-symmetric 
solitons with chemical potential fi — —10. The parameters 
are Vo = —4, go = — 1. Solid circles show the OL potential. 



chemical potential. Note that expression ((14]) is an even 
function of x, hence the soliton is symmetric. On the 
contrary, expression (|15|) is an odd function of x, which 
varies ~ x at x — > 0, yielding an antisymmetric soliton. 
In either case, the matter-wave densities are even func- 
tions of x. The exact soliton solution, cf> n , possess n — 1 
density nodes [see Figs. 2 and 3], and from Eq. (0) it 
can be concluded that the larger n, the larger the energy 
of the corresponding BEC state. Thus one may conclude 
that 4>i corresponds to the ground state, 4> n correspond- 
ing to the (n — l)-th excited states. By comparing the 
exact soliton solution 4>\ with the ground-state solution of 
the same GPE, obtained in a numerical form by means of 
the imaginary-time method, we find that <p\ is identical to 
the ground state when Vo < 0. However, <j>\ is not always 
the ground state when Vo > (for instance, <f>\ remains 
the ground-state solution at p < — 9 for Vq = 4). On the 
other hand, for the one-dimensional linear Schrddinger 
equation, it is well-known that localized states with dif- 
ferent energy eigenvalues are orthogonal. Here we find 
that the localized states of the nonlinear GPE are not 
orthogonal. 

From Figs. 1 and 2 it can be found that the widths 
of solitons are proportional to the widths of the respec- 
tive nonlinearity-modulation profiles, g{x). This is un- 
derstandable because both the widths of the solitons and 
g(x) profiles are determined by p(x), see Eqs. ^ and 
([7]). That is, the more rapidly p(x) approaches infin- 
ity, the narrower the solitons and g(x) distributions are. 
Further, it can be shown that the widths of the exact 
solitons are always larger than those of the respective 
modulation profiles. To analyze this point in a simple 
form, we here take the case of Vq — 0. In this case, 
p < p co = 0, p ~ exp(y/—px), so that g ~ exp(—6y/—px) 
and R ~ 1 — exp(— 1\J— px) at x > 0; thus the width of 
the soliton is about three times larger than that of the 
p(x) modulation. Since exact solitons in the left col- 
umn of Fig. 2 are confined mainly to a single OL cell, 
they can be called fundamental gap solitons [19| , whereas 
the right column displays broader gap solitons, alias gap 



FIG. 4: (Color online) Normalized exact asymmetric soli- 
tons (frnNn 1 ^ 2 , where N n is the corresponding norm (scaled 
number of atoms). Inset: the corresponding asymmetric 
nonlinearity-modulation profile. Here, we set go — — 1, with 
other parameters given in the text. Solid circles show the OL 
potential. 



waves [41] . 

Another noteworthy point is that, although the cutoff 
chemical potential, p co , is an even function of Vo, the 
nonlinearity-modulation function is not. From Eqs. ([5]) 
and (fT"3l) . it follows that the width of g(x) corresponding 
to Vq > is smaller than that for — Vq, see a typical 
comparison in Fig. 3. Thus, the effective nonlinearity 
in the case of — Vb is stronger than for Vb, at the same 
go- On the other hand, the effective potential is attractive 
(repulsive) for negative (positive) V> for fundamental gap 
solitons. For these reasons, the number of atoms in the 
exact fundamental gap solitons with Vo < is smaller 
than for Vo > 0, as shown by Figs. 2 and 3. 

Similarly, exact asymmetric solitons to Eq. ^) can 
be constructed if we let c\c 2 ^ in Eq. (fTTj) . For p 
to be sign-definite, we again need p < p co , and C\,C2 
should be carefully chosen. As a generic example, we 
take p = —5, Vb = 4, and C2 = go = — ci = — L In such 
a case, p(x) = MathieuC(— 5, 4, x) — MathieuS(— 5, 4, x), 
the asymmetric modulation profile is given by Eq. (J7J, 
and X(x) — J Q p~ 2 (s)ds. To meet the boundary condi- 
tions <j){x — > ±oo) = 0, constants A and X in Eq. (fill) 
should satisfy 



A[A(+oo) - X(-oo)] = 2nK(V2/2), 

X = XX(-oo)+K(V2/2), 



(16) 



where n = 1, 2, 3, • • • . 

The exact asymmetric solitons can be found when sub- 
stituting values (|T6| into Eqs. (|6]) and (|T0|) . The repre- 
sentative profiles of the solitons, together with the cor- 
responding asymmetric nonlinearity-modulation profile, 
are displayed in Fig. 4. Similar to the exact symmetric 
and antisymmetric solutions, different solitons with the 
same chemical potential are not orthogonal. 
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B. The case of E > 

The exact solution to Eq. ([5]) with E > is given 
by Eq. (jTDJ) . To construct the exact localized solutions, 
p(x) should approaches infinity as \x\ — > oo, so that the 
function X(x) is bounded; it can be shown that such a 
requirement may be realized when the chemical poten- 
tial p falls into the bandgaps of the spectrum induced by 
the OL potential of the corresponding linear Schrodinger 
equation. The exact solution for p is given by Eq. (|9|). 
Due to the nonzero value of the corresponding Wron- 
skian, (pi'tf2 — t P\ { P2 = — 1, ViO^) arL d tpzix) are not 
zero at the same position, and p is always sign-definite. 
However, when p does not belong to the semi-infinite 
bandgap, there exist several points where p is very close 
to zero, making the strength of the nonlinearity very large 
(this region is very narrow, and p looks like the delta 
function), which we do not consider here. We are rather 
interested in the case of p < p co . 

To meet the boundary condition, we need X[X(oo) — 
X(— oo)] = 2nK(m), n = 1, 2, 3, • • • . Because A > \f~E~, 
an inequality ensues from here, 

*(^)>™*. (17, 

From Eq. (IT7|) it follows that n > n max = [X(oo) — 
X{—oo)]\E/'K. Thus, unlike the case of E — 0, where 
n = 1,2,3,- ••, here the first several values of n may 
disappear. For example, if n max — 2.5, then actual val- 
ues which give rise to the solitons are n — 3,4,5, 
However, we find that, in the semi-infinite bandgap, 
i-max < 1, regardless of values of E, a, /?, and 7. That is 
to say, there is still an infinite number of exact solitons 
sharing the same chemical potential. The exact soliton 
solutions are given by Eqs. ©, ©, and (fTU)) . with the 
nonlinearity given by Eq. ([7]). 

C. The case of E < 

In this case, the sign-definite p exists when the chem- 
ical potential is in the semi-infinite bandgap, that is, 
p < p co , and the real constants a, /3, and 7 should be 
carefully chosen. To meet the boundary condition, we 
need A[X(oo) -X(-oo)] =2nK(m),n= 1,2,3, Un- 
like the case of E > 0, where inequality (|17p must be 
satisfied, here there is no restriction on n. That is to say, 
there is still an infinite number of exact solitons sharing 
the same chemical potential. 

D. Discussion 

Thus far, we have demonstrated above that an infinite 
number of exact soliton solutions can be constructed in 
the model with the OL potential, which share the same 
values of the chemical potential. These solutions exist 



in the semi-infinite bandgap, in accordance with the fact 
that families of gap solitons [n, for example, in Eqs. (|14p 
and (|15p , denotes the family's index] can be found in the 
semi-infinite bandgap when the attractive nonlinearity 
is spatially homogeneous [HI, [42[ • The same model also 
supports gap solitons in finite bandgaps; we are not going 
to discuss exact solitons in those bandgaps because the 
exact spatially modulated nonlinearity mimics the delta 
function, which (i) may be hard to realize in experiments, 
and (ii) the corresponding profile of exact solitons are 
irregular. 

We did not consider the repulsive nonlinearity here. 
The reason is that, for the repulsive nonlinearity, go > 
0, the nontrivial solution to Eq. ([5]) is U(x) = 
^J2{E~- A 2 )/. 9o sn(AX - X , ^E/\ 2 — 1), where A 2 < 
E < 2A 2 . To meet the boundary condition, we must 
demand A[X(oo) - X(-oo)] = 2nK{yjE/\ 2 - 1)), n = 
1,2,3, from which it follows that n < n max = 
[X(oo) — X(—oc)]"\/E/it. For the chemical potential 
falling into the semi- infinite bandgap, n max < 1< There- 
fore, there are no exact solitons in the semi-infinite 
bandgap for the spatially modulated repulsive nonlinear- 
ity, just like in the case of the spatially uniform repulsive 
nonlinearity [l|. 

IV. NUMERICALLY FOUND SOLITONS AND 
THEIR STABILITY 

In Sec. Ill, we were able to find only discrete sets of 
exact soliton solutions for the given nonlinearity. Here 
we consider more general matter-wave solitons with dif- 
ferent values of the chemical potential in the OL poten- 
tial, when the localized nonlinearity-modulation profile is 
fixed. That is, we aim to find solitons in the framework 
of equation 

^ = -^ + 2V co S (2^- MathieuC ^ Q) x)6 , (18) 

with po < p co and, generally speaking, p 7^ po, where p co 
is given by Eq. (fT2|) . We focus here only on symmetric 
soliton solutions in the semi-infinite bandgap, that is, 
p < p co . 

Exact solutions to Eq. (fT8|) have been found above for 
p = Pq. Thus, using the exact solutions as an initial 
guess, one can find more general solitons by means of 
the numerical relaxation method. There are two different 
cases, which we define as I and II, with g(x) localized, re- 
spectively, around a peak or bottom of the OL potential, 
for positive or negative Vo- In either case, solitons can be 
found in the semi-infinite bandgap, regardless of the value 
of ^o- In case II, the number of atoms is a monotonously 
decreasing function of p, just like in the case of the NLSE 
with the spatially uniform attractive nonlinearity. How- 
ever, the situation is quite different in case I. For the 
first family solitons, we find that the number of atoms 
at first decreases and then increases with the increase of 
p, see Fig. 5(b), whereas for other soliton families, the 
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FIG. 5: (Color online) (a) The first family of numerically 
found solitons. Solid (red), dashed (green), and dotted (blue) 
lines represent the solitons with chemical potentials fi = —10, 
li = —5, and fi = —4.5, respectively. The dashed-dotted line 
shows the nonlinearity-modulation profile, with /xo = — 10 and 
Vo — 4, see Eq. (|18[) (this corresponds to what is defined as 
case I in the text). Open circles show the OL potential, (b) 
The number of atoms versus the chemical potential, (c) and 
(d): The same as (a) and (b), except that the solitons are 
from the second family. 




FIG. 6: (Color online) The same as Fig. 5, except Vo = —4 
(which corresponds to case II, as defined in the text). 



atom number is a monotonously increasing function of 
fx, see Fig. 5(d). These types of the dependences have 
obvious implications for the solitons' stability, as per the 
Vakhitov-Kolokolov criterion [43[ , see below. For all soli- 
ton families in case I, when \i approaches the cutoff value, 
most atoms are located in wells of the OL potential adja- 
cent to the region where the nonlinearity is concentrated. 
In other words, the solitons are confined to one or two 
OL cells in case I, while in case II they are trapped in a 
single cell. 

Obviously, the stability of the solitons must be investi- 
gated too. To this end, we first employ the linear-stability 
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FIG. 7: (Color online) The largest instability growth rate 
versus the chemical potential. Solid lines, circles, and squares 
pertain to the first, second, and third families of the solitons, 
respectively. The parameters in (a) and (b) are the same as 
in Figs. 5 and Fig. 6, respectively. 




FIG. 8: (Color online) The evolution of an unstable soliton. 
The nonlinearity-modulation function and OL potential are 
the same as in Fig. 5(a). The initial condition is taken as per 
the numerically calculated solution [the dashed (green) line in 
Fig. 5(a)], mixed with a random (white-nose) perturbation. 



analysis. Substituting a perturbed solution, ijj(x,t} — 
[<fr(x) + u(x) exp(iSt) + v*(x) exp(— iS*t)] exp(—ifit), into 
Eq. |lj and linearizing it around the unperturbed one, 
4>{x), we arrive at the eigenvalue problem, 



C 



-C 



= 6 



(19) 



with operator C — d 2 /dx 2 +fi— 2Vb cos(2x) — 2g<p 2 . Here 
g = — MathieuC(^o, Vq, x)~ 6 . The soliton is unstable if 
any eigenvalue d has an imaginary part. 

Results of numerical calculations displayed in Fig. 7 
reveal that, in case II, the first and second families of soli- 
tons are stable against small perturbations, while higher- 
order soliton families are unstable. On the other hand, 
for case I, stable solitons emerge only in the first family, 
when the chemical potential is small enough, so that the 
solitons are very narrow, and the nonlinearity is strong 
enough to sustain solitons in the presence of the locally 
repulsive OL potential. Similar conclusions concerning 
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the stability of solitons supported by the competing (lo- 
cally attractive/repulsive) linear and nonlinear potentials 
were reported in Ref . [35[ . Comparing the date displayed 
in Fig. 7 with panels (b) and (d) in Figs. 5 and 6, we con- 
clude that the Vakhitov-Kolokolov criterion (dN/ dfi < 
as the necessary criterion for the stability of solitons sup- 
ported by the attractive nonlinearity [431 ]) holds in the 
present model. We have also checked the stability of 
exact soliton solutions. The results are in qualitatively 
agreement with those shown in Fig. 7. 

The stability of the solitons was further checked by 
direct numerical simulations of Eq. (TTJ) , producing results 
in agreement with the predictions of the linear-stability 
analysis. In particular, the solitons from the first family 
in the unstable region originally exhibit a quasi-stable 
evolution and then decay, with a larger part of the atom 
number located in a neighboring well of the OL, see Fig. 
8, while other unstable solitons quickly decay into noise. 

Although we have displayed here the results of the 
stability investigation only for two special nonlinearities, 
similar conclusions hold for other values of fj,o an d Vq as 
well. The asymmetric solitons too demonstrate a similar 
behavior. 



V. THE COMPOSITION RELATION BETWEEN 
SOLITONS AND NONLINEAR BLOCH WAVES 




FIG. 9: (Color online) The composition relation between 
NBWs and fundamental solitons. Opened circles in (a) and 
(b) represents fundamental solitons from the first family, while 
in (c) and (d) the circles denote fundamental solitons from the 
second family. Solid lines in (a) and (c) are NBWs found at 
the center of the Brillouin zone, while in (b) and (d) the solid 
lines depict the NBWs at the edge of Brillouin zone. The 
profiles of the numerically calculated NBWs completely over- 
lap with the expressions given by Eqs. fl 2 1 p and (122[) . The 
periodic nonlinearity-modulation profile is given by Eq. (|20|) . 
with go = — 1, fi — jito = —25, and Vo = —4. Solid circles 
show the OL potential. 



From Fig. 6, one can conclude find that the solitons 
and corresponding g(x) modulation profiles are confined 
to a single cell. Then, it may be interesting to form 
a spatially-periodic nonlinearity pattern, by placing the 
same local profiles of g(x) into other wells of the OL po- 
tential. In such a case, the system may admit not only 
the gap solitons, but also NBWs (nonlinear Bloch waves). 
For the NLSE with the spatially uniforms nonlinearity, 
the intuitive concept of the NBWs built as chains of fun- 
damental gap solitons has been recently justified in Ref. 
19], which has produced a composition relation between 
NBWs and fundamental solitons, although the relation 
cannot be expressed in a sufficiently simple mathemati- 
cal form. 

In this section, we demonstrate that the composition 
relation is also numerically valid in the GPE with the 
spatially periodic nonlinearity. To this end, we consider 
the following periodic nonlinearity-modulation pattern: 

9p{x) = ^ MathicuC(/^V , x - m7r)6 > ( 2 °) 

where m = 0, ±1, ±2, • • • , and the summation is per- 
formed over cells of the OL potential. 

For many values of /x, the single-peak modulation pro- 
file for g(x) given by Eq. (fT"3|) . and the respective soliton 
given by Eqs. (fl"4")l and (IT5*1) . can be confined to a single 
OL cell (here we focus on the symmetric case). For in- 
stance, g(x) and the soliton solution for fi = /Jo = —25 
at Vo = —4 meet this condition. In such cases, adja- 
cent solitons practically do not overlap, hence forces of 



the interaction between them in the periodic configura- 
tion are negligible. Therefore, one may try to represent 
NBWs as chains of fundamental solitons with identical 
or alternating signs: 

{<f>NBw) 1 = 224>n(x - rrnr), (21) 



{4>NBW) 2 = ^(-l) m <Mz - 7777T), (22) 

m 

where (j) n is given by Eq. (p3J or Eq. ([15]). The NBW 
described by Eq. (j2"Tj) is located at the center of the re- 
spective Brillouin zone, while that given by Eq. (|22|) is at 
its edge. The conjectured composition relation between 
the NBW and fundamental solitons was checked numer- 
ically for the first and second soliton families, as shown 
in Fig. 9. For other families of solitons, the composition 
relation also holds, in the same sense. 

When the soliton width exceeds the minimum period- 
icity of the OL, the composition relation between NBWs 
and fundamental solitons, as given by Eqs. (|2"Tj) and (|22l) . 
remains essentially valid, as long as the soliton widths 
are smaller than two OL periods, see a typical example 
in Fig. 10. In this case, in the cell where the soliton 
and NBW coexist, they overlap only in the central part 
of the solitons. The size of the region where the pro- 
files of the soliton and NBW overlap is determined by 
the soliton's width, hence it decreases with the increase 
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x 



FIG. 10: (Color online) The composition relation between 
NBWs and fundamental solitons. Dashed lines: the first fam- 
ily of the fundamental solitons. Solid lines: numerically ob- 
tained NBWs at the center of Brillouin zone. Opened cir- 
cles: NBWs given by Eq. (121 j) . The periodic nonlinearity- 
modulation profile, denoted by dash-dotted line in (a), is given 
by Eq. (f2T7fl with go = —1, Vb = 4, and fxo = —6. The chem- 
ical potentials of the solitons and NBWs in (a) and (b) are, 
respectively, /i = —6 and /i = —4.4. 

of the width. The composition relation effectively holds 
because the underlying GPE is almost linear in the re- 
gions where the two solitons overlaps. Numerical compu- 
tations also demonstrate that one can generate confined 
gap waves by putting several fundamental gap solitons 
together, with arbitrary combination of signs, cf. Refs. 
[R| E2] ■ However, as the soliton widths grows too large, 
the composition relation is no longer (numerically) valid, 
see Fig. 10(b). Thus, the numerically tested composite 
relation remains valid as long as the width of the in- 
dividual soliton does not exceed two OL periods (when 
the nonlinearity is spatially homogeneous, the composite 



relation remains valid as long as the width of the individ- 
ual soliton does not exceed one OL period, as conjectured 
and verified in another context in Ref. |19|). 

VI. CONCLUSIONS 

We have constructed an infinite number of exact soli- 
ton solutions, both symmetric and asymmetric, in the 
model of the BEC with the OL potential and specially 
devised profiles of the spatial modulation of the local 
attractive nonlinearity. The chemical potential of the ex- 
act solutions falls into the semi-infinite bandgap. These 
solitons may coexist, with different energies, at common 
values of the chemical potential. 

Based on the explicit solutions, we have also found 
generic solitons families in the numerical form, fixing the 
nonlincarity-modulation profile. The stability of the nu- 
merically found solitons has been checked by means of 
the linear-stability analysis, and also using direct simu- 
lations. 

Finally, we have discussed the composition relation be- 
tween nonlinear Bloch waves and the fundamental gap 
solitons. We have demonstrated numerically that the 
composition relation is virtually exact when widths of 
the solitons do not exceed the double period of the OL. 
When the width of solitons exceeds the OL period, the 
size of the spatial region where the solitons and nonlin- 
ear Bloch waves overlap decreases with the increase of 
the soliton's width. 
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